function fig5_5b

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)

nx=500;

aa=1; ba= 1; ep=1/25;
t0 = sqrt(3) - log(2 + sqrt(3))/2;
a0 = cos(t0/ep - pi/4);
a1 = cos(t0/ep + pi/4);
aR = (aa*3^(1/4) - ba*a1*exp(-pi/(4*ep))) / (2*a0 - a1*exp(-pi/(2*ep)));
bR = exp(-pi/(4*ep))*(ba - aR*exp(-pi/(4*ep)));
A = 2^(5/6)*sqrt(pi)*aR;
B = bR*sqrt(pi)/2^(1/6);

xxL=linspace(-1,0.001,nx);
for ix=1:nx
    t = ( 1 - xxL(ix) )*sqrt( xxL(ix)^2 - 2*xxL(ix) ) /2 - log( 1 - xxL(ix) + sqrt( xxL(ix)^2 - 2*xxL(ix) ) )/2;
    qa = xxL(ix)*(2-xxL(ix));
    yyL(ix) = (2*aR*cos(t/ep - pi/4) + bR*cos(t/ep + pi/4))/abs(qa)^(1/4);
end

xxR=linspace(0.001,1,nx);
for ix=1:nx
    k = (xxR(ix)-1)*sqrt(2*xxR(ix)-xxR(ix)^2)/2 + asin(xxR(ix)-1)/2 + pi/4;
    qa = xxR(ix)*(2-xxR(ix));
    yyR(ix) = (aR*exp(-k/ep) + bR*exp(k/ep))/qa^(1/4);
end


%  y'' + p(x)y' + q(x)y= f(x)   for xL < x < xR
% set boundary conditions
xL=-1; yL=1;
xR=1; yR=1;

x=linspace(-1,1,nx);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
a=zeros(1,nx-2); b=zeros(1,nx-2); c=zeros(1,nx-2);
for i=1:nx-2
    a(i)=-2+h*h*q(x(i+1),ep);
    b(i)=1-0.5*h*p(x(i+1),ep);
    c(i)=1+0.5*h*p(x(i+1),ep);
    f(i)=h*h*rhs(x(i+1),ep);
end
f(1)=f(1)-yL*b(1);
f(nx-2)=f(nx-2)-yR*c(nx-2);
% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);
y=[yL, y, yR];

hold on
box on
grid on

plot(x,y,'-r','Linewidth',1.1)

plot(xxL,yyL,'--b','Linewidth',1.1)

plot(xxR,yyR,'--k','Linewidth',1.1)


axis([-1 1 -3 6])
set(gca,'ytick',[-2 0 2 4 6])
set(gca,'YTickLabel',{'$-2$';'$0$';'2';'4';'6'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Solution','Interpreter','latex','FontSize',14)
text(-0.95,5,'(b)','Interpreter','latex','FontSize',14)
legend(' $\, y$',' $\, y_L$',' $\, y_R$','Location',[0.727 0.623 0.115 0.146],'FontSize',13,'Interpreter','latex')

%exportgraphics(gcf,'/Users/mark/Desktop/fig44b.pdf')


function g=q(x,ep)
g=-x*(2-x)/ep^2;

function g=p(x,ep)
g=0;

function g=rhs(x,ep)
g=0;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
    y(j) = y(j) - v(j)*y(j+1);
end


